#include <stdio.h>
#include <stdlib.h>
#include "gmp.h"
#include <pthread.h>

#define ITER            25
#define PRECISION       95000000
#define BASE            10
#define DIGITS          30000000
#define N_THREADS       3

int iter;

/* Main Variables */
mpf_t a[ITER + 1], b[ITER + 1], p[ITER + 1];

/* MUTEX */
pthread_mutex_t amutex[ITER + 1];
pthread_mutex_t bmutex[ITER + 1];
pthread_mutex_t pmutex[ITER + 1];

/* Headers of the functions that calculates a, b and p variables */
void *calc_a();
void *calc_b();
void *calc_p();

int
main(int argc, char *argv[])
{
    pthread_t threads[N_THREADS];
    long i;

    fprintf(stderr, "Calculando...\n");

    if(argc == 2)
        iter = iter < 25 ? atoi(argv[1]) : ITER;
    else iter = ITER;

    /* Initializes all of the mutexes */
    for(i = 0; i <= iter; i++)
    {
        pthread_mutex_init(&amutex[i], NULL);
        pthread_mutex_init(&bmutex[i], NULL);
        pthread_mutex_init(&pmutex[i], NULL);
    }

    /* Defines the precision and initializes all the variables */
    mpf_set_default_prec(PRECISION);
    for(i = 0; i <= iter; i++)
        mpf_inits(a[i], b[i], p[i], NULL);

    /* a = sqrt(2) */
    mpf_init_set_ui(a[0], 2);
    mpf_sqrt(a[0], a[0]);

    /* b = 0 */
    mpf_init_set_ui(b[0], 0);

    /* p = 2 + sqrt(2) */
    mpf_init_set_ui(p[0], 2);
    mpf_add(p[0], p[0], a[0]);

    /* Lock all the mutexes except those of the initial variables */
    for(i = 1; i <= iter; i++)
    {
        pthread_mutex_lock(&amutex[i]);
        pthread_mutex_lock(&bmutex[i]);
        pthread_mutex_lock(&pmutex[i]);
    }

    /* Creates the threads and set each to one function */
    pthread_create(&threads[0], NULL, calc_a, NULL);
    pthread_create(&threads[1], NULL, calc_b, NULL);
    pthread_create(&threads[2], NULL, calc_p, NULL);

    /* Waits for all of the threads to finish */
    for(i = 0; i < N_THREADS; i++)
        pthread_join(threads[i], NULL);

    /* Destroy all of the mutexes freeing memory */
    for(i = 0; i <= iter; i++)
    {
        pthread_mutex_destroy(&amutex[i]);
        pthread_mutex_destroy(&bmutex[i]);
        pthread_mutex_destroy(&pmutex[i]);
    }

    /* Prints the final result */
    mpf_out_str(stdout, BASE, DIGITS, p[iter]);

    /* Destroy all of the variables freeing memory */
    for(i = 0; i <= iter; i++)
        mpf_clears(a[i], b[i], p[i], NULL);

    return 0;
}

void *
calc_a()
{
    /* Inicialization of the variables used in this function */
    mpf_t sqrt_a, inv_sqrt, sum1;
    mpf_set_default_prec(PRECISION);
    mpf_inits(sqrt_a, inv_sqrt, sum1, NULL);

    int i;

    /* an1 = (sqrt(an) + (1 / sqrt(an))) / 2; */
    for(i = 0; i < iter; i++)
    {
        mpf_sqrt(sqrt_a, a[i]);
        mpf_ui_div(inv_sqrt, 1, sqrt_a);
        mpf_add(sum1, sqrt_a, inv_sqrt);
        mpf_div_ui(a[i + 1], sum1, 2);

        /* Unlocks the variable meaning that a[i+1] is ready for use*/
        pthread_mutex_unlock(&amutex[i + 1]);
    }

    mpf_clears(sqrt_a, inv_sqrt, sum1, NULL);
    return;
}

void *
calc_b()
{
    /* Inicialization of the variables used in this function */
    mpf_t sum1, sqrt_a, mult, sum2;
    mpf_set_default_prec(PRECISION);
    mpf_inits(sum1, sqrt_a, mult, sum2, NULL);

    int i;

    /* bn1 = ((1 + bn) * sqrt(an)) / (an + bn) */
    for(i = 0; i < iter; i++)
    {
        /* Test if the variable a[i] is ready */
        pthread_mutex_lock(&amutex[i]);
        pthread_mutex_unlock(&amutex[i]);

        mpf_add_ui(sum1, b[i], 1);
        mpf_sqrt(sqrt_a, a[i]);
        mpf_mul(mult, sum1, sqrt_a);
        mpf_add(sum2, a[i], b[i]);
        mpf_div(b[i + 1], mult, sum2);

        /* Unlocks the variable meaning that b[i+1] is ready for use*/
        pthread_mutex_unlock(&bmutex[i + 1]);
    }

    mpf_clears(sum1, sqrt_a, mult, sum2, NULL);
    return;
}

void *
calc_p()
{
    /* Inicialization of the variables used in this function */
    mpf_t sum1, res1, res2, sum2;
    mpf_set_default_prec(PRECISION);
    mpf_inits(sum1, res1, res2, sum2, NULL);

    int i;

    /* p = ((1 + a) * p * b) / (1 + b) */
    for(i = 0; i < iter; i++)
    {
        /* Test if the variable a[i+1] and b[i+1] is ready */
        pthread_mutex_lock(&amutex[i + 1]);
        pthread_mutex_lock(&bmutex[i + 1]);
        pthread_mutex_unlock(&bmutex[i + 1]);
        pthread_mutex_unlock(&amutex[i + 1]);

        mpf_add_ui(sum1, a[i + 1], 1);
        mpf_mul(res1, sum1, p[i]);
        mpf_mul(res2, res1, b[i + 1]);
        mpf_add_ui(sum2, b[i + 1], 1);
        mpf_div(p[i + 1], res2, sum2);

        /* Unlocks the variable meaning that p[i+1] is ready for use*/
        pthread_mutex_unlock(&pmutex[i + 1]);
    }

    mpf_clears(sum1, res1, res2, sum2, NULL);
    return;
}
